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Motivated by the initial-boundary value problem for the Einstein equations, we propose a defini- 
tion of symmetric hyperbohcity for systems of evolution equations that are first order in time but 
second order in space. This can be used to impose constraint-preserving boundary conditions. The 
general methods are illustrated in detail in the toy model of electromagnetism. 



I. INTRODUCTION 

In attempting to solve the Einstein equations numeri- 
cally as an initial or initial-boundary value problem, it is 
important to start with a well-posed continuum problem. 
For systems of first-order equations, there are two useful 
criteria for well-posedness: Strong hyperbolicity is a nec- 
essary and sufficient criterion for the initial value prob- 
lem (without boundaries or with periodic boundaries) to 
be well-posed. Roughly speaking, it means that the sys- 
tem possesses a complete set of characteristic variables, 
which propagate with definite speeds. Symmetric hyper- 
bolicity implies strong hyperbolicity, and can be used to 
formulate a well-posed initial-boundary value problem. 
Roughly speaking, it means that the principal part of the 
system possesses an (unphysical) energy whose growth 
due to boundary terms can be expressed as "incoming 
mode energy minus outgoing mode energy" . 

Much less is known about how to make systems of 
second-order equations well-posed. Here we concentrate 
on systems which have been brought into first order in 
time but remain second order in space. The simplest 
example of such a system is the wave equation, written 
as 



M = n, 



(1) 
(2) 



The Arnowitt-Deser-Misner (ADM) formulation of the 
Einstein equations has a similar form, with the 3-metric 
corresponding to and the extrinsic curvature to H. 

In this paper we discuss three methods for investigating 
and proving well-posedness for such systems: we review 
two known methods in order to put our work into context, 
and then propose one which we believe to be new and 
advantageous. 

The first old method, (differential) reduction to first 
order, introduces the spatial derivatives as auxiliary vari- 
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ables. For the wave equation above, these could be called 
di = di4>. One then has a larger but first-order system, 
which one can hope to make strongly or symmetric hy- 
perbolic. A drawback of this method is that its solution 
space is larger than that of the original system: it is 
equivalent to the original system only if the three auxil- 
iary constraints Ci = di — dicf) vanish. Therefore proving 
strong or symmetric hyperbolicity for the first-order re- 
duction does not prove well-posedness for the original 
second-order system except in simple cases where the 
auxiliary constraints evolve trivially. (These include the 
wave equation and electromagnetism, but not the full 
Einstein equations.) 

A second method called pseudo-differential reduction 
to first order Fourier-transforms the PDE system in 
space. If it is linear with constant coefficients, one effec- 
tively obtains a system of decoupled ODEs, one for each 
wave number LOi. One can then show a version of strong 
hyperbolicity for this system. The key point is that the 
pseudo-differential reduction does not introduce any aux- 
iliary variables. Proving strong hyperbolicity therefore 
also proves well-posedness of the initial value problem 
for the original second-order system 0. This method 
has been applied to a formulation of the Einstein equa- 
tions in . A drawback is that as the Fourier transform 
is nonlocal, one cannot find a locally conserved energy 
in this way, and so cannot address the initial-boundary 
value problem. 

In order to avoid the drawbacks of these two ap- 
proaches, we propose a fairly obvious generalisation of 
the concept of symmetric hyperbolicity to systems that 
are second order in space and first order in time. In the 
example of the wave equation, we define characteristic 
variables in a given direction Ui as n±rt*9j0, and an en- 
ergy density 11^ -|- {di4>)'^. The energy and its associated 
flux can be written in terms of characteristic variables, 
and so we can impose maximally dissipative boundary 
conditions similar to those for a first-order symmetric 
hyperbolic system. 

The paper is organised as follows: In Sec. we re- 
view strong and symmetric hyperbolicity, and the initial- 
boundary value problem, for first-order systems. In 
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Sec. mil we review pseudo-differential reduction to first 
order. In Sec. IIVI we present our concept of symmetric 
fiyperbolicity for a system that is second order in space 
and first order in time. In Sec. we review the initial- 
boundary value problem for systems with differential con- 
straints, such as the Maxwell or Einstein equations. This 
discussion applies both to first and to second order sym- 
metric hyperbolic systems. 

Finally, in Sec. lVII we compare the three different meth- 
ods by applying them to the Maxwell equations: re- 
duction to first order with auxiliary variables, pseudo- 
differential reduction to first order, and our method of 
finding characteristic variables and an energy directly 
for the second-order system. We formulate constraint- 
preserving boundary conditions, show that the initial- 
boundary value problem for the second-order Maxwell 
equations can be made well-posed for "Neumann" and 
"Dirichlet" boundary conditions, and conjecture that it 
is well-posed for a one-parameter family of boundary con- 
ditions which interpolates between those two via "Som- 
merfeld" boundary conditions. 

II. FIRST-ORDER SYSTEMS 

In this section we summarise some general results 
on the well-posedness of systems of first-order evolution 
equations that are proved for example in We begin 
by considering linear systems with constant coefficients, 
and at the end of the section appeal to standard methods 
for generalising these results to quasilinear systems. (For 
other reviews, see 

A. Strong hyperbolicity 

Consider the system of first-order linear evolution 
equations with constant coefficients 

dtu = P'diU + Qu. (3) 

Here u is a vector of variables, P* and Q are square matri- 
ces, and i ranges over the spatial directions. The initial- 
value problem for such a system is called well-posed if 
a solution exists, is unique and depends continuously on 
the initial data u(x, 0). This is equivalent to there being 
an estimate 

|K.,t)||<i^e"*|K.,0)|| (4) 

with respect to an L2 norm, where the constants K and 
a do not depend on the initial data. To simplify the 
presentation, in the following we shall assume that Q = 0. 
Leaving Q in the calculation would show that it influences 
the value of a, but not the existence of the estimate itself. 

The main result is that the existence of an estimate 
Q is equivalent to the system being strongly hyperbolic: 
for any unit covector n^, the matrix P„ = riiP^ has only 
real eigenvalues and a complete set of real eigenvectors. 



The idea of the proof of this theorem is to go to Fourier 
space. Let u{LU,t) be the Fourier transform of u(x,t): 

ui.,t)^jji.,t)e-'-^^. (5) 

Parseval's theorem states that \\u\\x — ^^'^ the 

estimate (0} is equivalent to 

\\u{-,t)\\<Ke"'M-M\ (6) 

The system in x and t becomes an ODE system for each 
value of uji : 

dtu = i|w|P„M, (7) 

where now = Wi/lajj. Now, if P„ has complex eigen- 
values, then a depends on the largest present in the 
initial data, and hence on the initial data. If P„ lacks a 
complete set of eigenvectors, the normal form of P„ con- 
tains at least one Jordan block of size k > 1. One can 
show that then \u{uj,t)\ is bounded only by |a;t|'^~-^, so 
that now K depends on the largest value of |a;| present 
in the initial data. 

On the other hand, if Pn is diagonalisable with real 
eigenvalues, let A„ be the diagonal matrix of eigenvalues, 
and T„ be a matrix of the corresponding right (column) 
eigenvectors of P„, 

PnTn = TnAn- (8) 

(Note that r„ is not unique.) The Fourier characteristic 
variables 

U{co,t) = T-'uiLo,t) (9) 

obey 

dtU = t\uj\A,,U. (10) 

This gives 

\U{Lu,t)\^\U{u;,0)\, (11) 

and so 

\u{u;,t)\<\T^\\T-^\\u{u;,0)\. (12) 

If Tn and T^^ are continuous in rii, we obtain © and 
hence (gj, with a = and K = sup(|T„| \T-^\). 

We now go back to real space, and define the vector 
U of characteristic variables with respect to an arbitrary 
fixed direction rii by 

U = T-\. (13) 

The characteristic variables always depend on the direc- 
tion Ui, but following standard notation we do not in- 
dicate this explicitly. We define the derivative in the 
direction rii as 

dn = n'9„ (14) 
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and the projector ^ into the space normal to rii by 

q'^=5\j-n'nj. (15) 

(Note that we need a spatial metric to define n*.) We 
also define the shorthand index notation A for a tensor 
index i that has been projected with g*^, for example in 

p^^d, = p^dA (16) 

With this notation we can write 

P'd^^Pndn+P^dA. (17) 

The characteristic variables obey 

dtU = AndnU + {T,-;^P^Tn)dAU. (18) 

Strong hyperbolicity implies the existence of a sym- 
metriser, which is a Hermitian matrix iJ„ for every unit 
covector rii that obeys 

(ff„P„)t = H„P„ (19) 

The general matrix Hn with this property can be con- 
structed as 

F„ = (T-i)ti?„r-i, (20) 

where Bn is Hermitian, positive definite, and commutes 
with A„. The general matrix B„ with this last property 
is block-diagonal, with each block corresponding to one 
(possibly multiple) eigenvalue in A„. For given T„, Hn 
is parameterised by -B„. Conversely, it can be parame- 
terised by T„ with Bn = I. If the field equations do not 
single out a preferred direction n^, then A„ do not de- 
pend on Hi, and it is natural to choose r„ so that Bn is 
also independent of rii. We then write B and A instead 
of Bn and A„. 



B. Symmetric hyperbolicity 

We shall now see how well-posedness can be main- 
tained in the presence of boundaries in space. With- 
out loss of generality we assume that these boundaries 
are fixed (at constant values of the coordinates a;'.) If 
a symmetriser Hn = H exists that is independent of , 
the system is called symmetric hyperbolic. In this case 
the system admits a positive definite energy 

I edV, e = v)Hu. (21) 
Jn 

This means that E is conserved up to boundary terms: 
Assuming still that H and are constant, and with 
HP' = {HP')\ we have 

dte = d,{u^HP'u), (22) 



and so by Gauss' law 
dE f 

— = / F'dS,, F' = v)HP'u, (23) 

where dSi is the surface element on the boundary dfl. 
This can be written as 

dE r 

— = / F" dS, P" = U'fBAU, (24) 
dt Jdn 

where is the normal to dfl. 

The boundary flux can be written schematically as 
"[/^ — C/^", or more precisely as 

P" = C/|P+A+C/+ + [/lP_A_t/_, (25) 

where [/+ are modes that come in across the boundary 
and U- are outgoing, and where B± and A-j- are the 
blocks corresponding to U±. Therefore B+A+ is positive 
definite and _B_A_ is negative definite (the meaning of 
the minus sign in "[/^ — f/^"). There may also be modes 
Uq whose propagation speed across the boundary is zero. 
In this case the boundary is called characteristic. These 
modes do not contribute to F". It is clear that boundary 
conditions must be imposed on the incoming variables 
[/_!_, and that no boundary conditions can be imposed 
on the variables U- and Uq. In the following we assume 
that the number of variables that propagate along the 
boundary is constant. 

Estimates of the type can be obtained from an 
energy E whose growth is suitably bounded. One can 
achieve this by imposing a maximally dissipative bound- 
ary condition 

U+ = LU- + /, (26) 

where / is given and the matrix L is sufficiently small. 
For f — this condition means that more energy goes 
out than comes in at every point on the boundary (not 
just integrated over the boundary). For / ^ given, 
the growth of E can be bounded from /. The required 
smallness condition on the matrix L is that the matrix 

L^B+K+L + B^K^ (27) 

is negative definite. For maximally dissipative bound- 
ary conditions to result in a well-posed initial-boundary 
value problem, the system has to be symmetric hyper- 
bolic. There are examples where imposing them on a 
strongly but not symmetric hyperbolic system leads to 
an ill-posed problem |^]. 

C. Variable coefficient and quasilinear systems 

The results on linear systems with constant coefficients 
can be generalised to linear systems with variable coeffi- 
cients 13 • It is sufficient that the system is strongly or 
symmetric hyperbolic at every point in space and time 
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(frozen coefncients approximation), and bounds can be 
obtained that are uniform over all space, all directions 
Hi, and a finite time interval. In obtaining an energy 
for the system with variable coefhcients, the time and 
space derivatives of the non-constant cocfhcients give 
rise to non-principal terms, but these do not affect well- 
posedness. 

The results for linear problems with variable coeffi- 
cients can be extended to nonlinear problems. In a first 
step, linearising a non-linear problem around a back- 
ground solution gives rise to a linear problem with vari- 
able coefficients. Clearly a necessary condition for the 
non-linear problem to be well-posed is that the linearised 
problem is well-posed. Furthermore, if the linearised 
problem is well-posed for all backgrounds, then it may be 
possible to prove well-posedness for the nonlinear prob- 
lem 0. 

Physically, the frozen coefficients, principal part only 
approximation means that we focus on possible instabili- 
ties in small high-frequency perturbations. A problem is 
ill-posed in practice because it allows perturbation modes 
whose growth rate is not bounded - the higher the fre- 
quency in space, the faster the growth in time. In nu- 
merical simulations such modes show up as instabilities 
at the grid frequency which grow faster at higher res- 
olution, thus destroying convergence at sufficiently late 
times and high resolutions. Well-posedness does not rule 
out instabilities (in the loose sense of rapidly growing 
perturbations, possibly violating the constraints) which 
are generated by non-principal terms and whose growth 
is bounded independently of spatial frequency. 



The Fourier transform of this system is 

dt$ = n, (28) 
dttl = -\u;\^. (29) 

Replacing (f) by ip = i\uj\(j) 0, we obtain the first-order 
pseudo-differential system 

dtip = i\uj\U, (30) 
dtfl = i\uj\i^. (31) 

Note that this is not the Fourier transform of a differen- 
tial system as P„ is independent of n^. This system is 
strongly hyperbolic, with Fourier characteristic variables 

U±=tl±ip. (32) 

which obey 

dtU± = ±i\uj\U±. (33) 

IV. SECOND ORDER IN SPACE, FIRST 
ORDER IN TIME 

We shall now show how the two key ideas for first-order 
systems - strong and symmetric hyperbolicity - can be 
adapted to a system that is first order in time but second 
order in space. 

A. Strong hyperbolicity 



III. PSEUDO-DIFFERENTIAL REDUCTION TO 
FIRST ORDER 



Pseudo-differential reduction to first order has recently 
been suggested as a technique for proving that certain 
second-order evolution systems are strongly hyperbolic 
0. We review it here only to see how it relates to the 
two other ways of dealing with a second-order system: 
reduction to a first-order system by introducing auxil- 
iary variables, and the second-order energy method we 
propose in this paper. 

A system of first order pseudo-differential equations 
is one whose Fourier transform in space is of the form 
(|7|) but where P„ is not given by the linear expression 
Pn = riiP"^ for a constant P*. The definition of strong hy- 
perbolicity of a system of evolution equations generalises 
to the pseudo-differential case, and strong hyperbolicity 
is again equivalent to well-posedness (in the absence of 
boundaries), because the proof using the ODE system 
O refers only to the matrix P„ but not directly to the 
vector of matrices P* . 

First-order pseudo-differential systems can be derived 
from a class of second-order differential systems. The 
simplest example is the wave equation, in the form H1I2|I . 



In analogy with H18|l we define a linear combination 
U oi u = {u, diu) to be a characteristic variable of the 
second order system with speed A in the direction ni if it 
obeys 



dtU — XdnU + transversal derivatives. 



(34) 



We illustrate this again for the scalar field model. The 
second order characteristic variables are 



U± = U±dn^ 
Ua = 



and obey 



dtU± = ±a„c/± + 9^f/A, 

dtUA - dA\{U+ + U-). 



(35) 
(36) 



(37) 
(38) 



Note that one can obtain the non-zero speed charac- 
teristic variables H35|l of the second order PDE system 
from those of the pseudo-differential system, (|32|l . be- 
cause both are eigenvectors of essentially the same matrix 
P„ . In the second-order system we neglect the transverse 
derivatives to find P„, while in the pseudo-differential 
system is defined as Wi/jwl, and so the transversal 
derivatives vanish by definition. 
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The definition 1)34(1 of characteristic variables for a 
second-order in space system has two consequences which 
are not apparent in the simple scalar field example. The 
first is that a transversal derivative Oau is by definition 
a zero-speed characteristic variable because we can inter- 
change derivatives: 

dtidAu) = dA{dtu). (39) 

The second is that all characteristic variables of a second- 
order system are unique only up to addition of transversal 
derivatives. This non-uniqueness will later be used up to 
write e and in terms of characteristic variables. 



B. Symmetric hyperbolicity 

We would like to generalise the energy to a sys- 
tem that is second order in space. As in the first order 
system, we require an energy density and flux that obey 
the conservation law 

dte = d,F\ (40) 

and can be expressed in characteristic variables. 

We demonstrate the second-order energy and maxi- 
mally dissipative boundary conditions for the scalar field. 
The scalar field admits a positive definite conserved en- 
ergy that is unique up to an overall scale, with energy 
density 

e - + (a,</>)2 (41) 

and fiux 

F' = 2n5>. (42) 
The energy can be written as 

e=\{Ul + Ul) + UAU^ (43) 

and the fiux in the direction Ui normal to the boundary 
can be written as 

= i (C/2 _ ul) . (44) 

This means that 

U+ = KU-+f (45) 

with |k| < 1 is a maximally dissipative boundary condi- 
tion for the system 1)112(1 . As the characteristic variables 
of the second-order system contain derivatives, this is not 
an algebraic but a first-order differential boundary con- 
dition, namely 

(1 - K)n + (1 + K)dn(j} - /. (46) 

We see that k = —1 corresponds to an (inhomogeneous) 
Dirichlet boundary condition on 11 and hence on cf), k = 1 



to Neumann, and k = to specifying the ingoing radia- 
tion (Sommerfeld). 

The conserved energy with e given by (|^ can be used 
to bound the L2 norms of 11 and dicf), but not ||0|| itself. 
To obtain a bound on \ \ as well, we can use the energy 
given by 

with a > 0. Clearly this gives us 

||n||<^i/2, \\d,cl>\\<E^/^ M\<a-'E'/^ (48) 

but with the disadvantage that E is not strictly con- 
served. Rather 

^ = / 2Ud'(j)dS, + 2a^ [ UcjjdV < 2aE (49) 
Jan Jn 

where we have assumed maximally dissipative boundary 
conditions to make the boundary term non-positive. This 
is an inefficient estimate, but a can be made arbitrarily 
small. Non-principal terms in the equations can be dealt 
with in a similar manner: they give rise to volume terms 
in dE/dt that can be estimated as a multiple of E itself. 
(Note that the issues of bounding \ \ and of dealing with 
non-principal terms are dealt with identically in second- 
order and first-order systems.) 



V. CONSTRAINT-PRESERVING BOUNDARY 
CONDITIONS 

Here we summarise a general method of dealing with 
a system of evolution equations that are subject to dif- 
ferential constraints (0 and references there). As the 
constraints are by assumption compatible with the evo- 
lution equations, they form a closed evolution system. 
We now assume that this constraint system and the main 
system are symmetric hyperbolic. The following discus- 
sion applies then equally to first-order and second-order 
systems. 

Rather than defining a notation for the general case, we 
consider the case where the non-zero speed characteristic 
variables come in pairs with speeds ±A, with A > 0, plus 
a number of zero-speed variables. (This is necessarily 
true if there is no preferred direction in space.) Focus 
on a pair of characteristic constraint variables C± which 
obey 

dtC± = ±XdnC± + . . . , (50) 

where . . . stands for transversal derivatives and any non- 
principal terms. As the constraints are compatible with 
the evolution equations, there exists a pair of characteris- 
tic variables of the main system U± with the same speeds, 
that is 

dtU± = ±\dnU± + . . . , (51) 
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such that 



be written as the evolution equations 



(52) 



The constraint energy is initially zero and we want it to 
remain zero. One way of enforcing that is to formally im- 
pose homogeneous maximally dissipative boundary con- 
ditions on the constraint system, that is 



C+ = kC- 
for \k\ < 1. If we now define 



(53) 



(54) 



the boundary condition on the constraint system 1)53(1 is 
equivalent to the boundary condition 



dnX 



(55) 



on the main system, taking into account the relation (|52|l . 
We now use the evolution equations for U± to replace 
dn with dt, and find an evolution equation for X which 
contains only transversal derivatives (derivatives parallel 
to the boundary): 



dtX 



(56) 



The trick is to interpret this as an evolution equation 
for an auxiliary variable X which is defined only on the 
boundary. (In general, there will be more than one pair of 
constraints, and a system of variables X.) The boundary 
condition 



kU^ +X 



(57) 



has the form of a maximally dissipative boundary condi- 
tion on U, but it is one properly speaking only if X can 
be treated as given. This is the case only if the boundary 
system (|56|l decouples from the bulk system. 

One can show that, with maximally dissipative bound- 
ary conditions, 



^ll^llao < K.WfWon + K2\\u\\an, 
at 

^A\u\\n < K^WXWon. 
at 



(58) 
(59) 



If the boundary system decouples from the bulk system, 
K2 — and the growth of the boundary solution and 
hence of the bulk solution is bounded by /. But a bound 
on ||w||n does not imply one on ||M||an, and so for K2 ^ 
well-posedness cannot be proved using these energy tech- 
niques. 



VI. THE MAXWELL EQUATIONS 

A. Field equations 

The Maxwell (electromagnetism, EM) equations in the 
presence of a charge density p and current density ji can 



dtE, = ~djd' A + d^d' Aj - Airj, 
subject to a constraint 

Ce = d'Ei - 4ttp = 0. 



(60) 
(61) 



(62) 



Here ^, ji and p can be treated as given functions, subject 
to the charge conservation law dtp + d^ji = 0. All indices 
are moved with the flat Cartesian spatial metric 6ij , and 
staggered double indices are summed over. 

The EM system would be a wave equation for Ai were 
it not for the grad div term did^Aj. This term can 
be eliminated by making the offending divergence a new 
variable. This property makes the EM system a nice toy 
model for the ADM formulation of the Einstein equations 
;8j. We therefore introduce F = d^Ai. This gives rise to 
the new constraint 



Cr = r - d'A. = 0. 



(63) 



We obtain an evolution equation for T by taking the di- 
vergence of the evolution equation for Ai, and add to it 
bCs with a free coefficient b. We also add adiCy to the 
evolution equation for Ei with a free coefficient a. The 
resulting system is 

dtA, = -E,-d,i>, (64) 
dtEi = -djd^Ai + adiT + (1 - a)d^d^ Aj - Attj,, (65) 
dtT = {h-l)diE' -hAnp-did'ip. (66) 

Note that for a = 0, F decouples and we recover the 
original system. For a — 1, Ai obeys a wave equation. 
The constraint system is 



dtCE 
dtCr 



ad.d'Cr, 
bC'E. 



(67) 
(68) 



In the following, we concentrate on the principal part. 
In the EM system, this means neglecting the terms con- 
taining p, ji and ip. (For the EM system the principal 
part approximation happens to be exact: it corresponds 
to vacuum electrodynamics in the gauge ^ = 0.) 

We now use the Maxwell equations to illustrate three 
methods of investigating well-posedness: differential re- 
duction to first order, pseudo-differential reduction, and 
making the original second-order system symmetric hy- 
perbolic. This is only for illustration: our preferred 
method is the last one. We have introduced the param- 
eters a and b to clarify the presentation, and because 
precisely analogous parameters exist for the ADM for- 
mulation of the Einstein equations. We shall see that a 
natural choice of these parameters would be a = 6 = 1. 
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Differential reduction to first order 



Field equations 



To reduce the EM main system to first order, we in- 
troduce the auxiUary variables dij = diAj. Because Ai 
docs not appear nndifFcrcntiatcd in the field equations, 
this definition itself does not appear as a new constraint, 
but its integrabihty condition does: 



0. 



(69) 



The only term in the second-order system that does not 
translate unambiguously into the new variables is did^Aj, 
which can be written either as did/ or d^dij. Wc there- 
fore paramcterise this choice with a parameter c. From 
the point of view of the new system, this can be described 
as adding the constraint —cCij-' to the evolution equation 
for Ei with a free coefficient c. The resulting system is 
(now neglecting non-principal terms as discussed above) 



dtdiA = —diEi 



dtEi 
dtT 



(70) 



-ff'dji + adiT + {l-a- c)did/ + cd^ij, (71) 
{b-l)diE\ (72) 



To reduce the constraint system to first order, we 
introduce the new variables C, = diCr, where now 
Cr = T — dj-' . As Cr does not appear in undifferen- 
tiated form, this definition itself does not appear as a 
constraint, but it gives rise to the integrabihty condition 



Cii 



diCj 



djCi = 0. 



The first-order constraint system is 

dtCE = ad'Ci + {l-c)d'Ci/, 

diC, = bd.CE, 

dtCijk = 0, 
dtCii = 0. 



2. Strong hyperbolicity 

We write the first-order system as 

dtu = PndnU + P"^d„Au. 



(73) 



(74) 
(75) 
(76) 
(77) 



(78) 



To diagonalisc P„ . we split dij into the tensor compo- 
nents dnn, dqq, dAn, dnA and dAB with respcct to a fixed 
direction rij, where the index n means a contraction with 
n*, the index pair qq means a contraction with g*^, and 
indices A, B, etc are the same as i, j etc, but denote a 
tensor that is transversal and tracefree on any index pair. 
(Because it is transversal, it is tracefree with respect to 
both Sij and qij.) Note that di^ = dnn + dqq. A similar 
tensor decomposition is used for the other variables and 
throughout this paper. We shall loosely refer to these ob- 
jects as scalars if they have no free (transversal) index. 



vectors if they have one, and tensors if they have more 
than one. 

By applying the same tensor decomposition to the 
evolution equations, we find P„ in this basis: it is 
block-diagonal with blocks corresponding to the scalars 
{dnn, dqq, En, T), transversc vectors {d„A,dAn,EA) and 
transverse tracefree tensors (d^s)- The free transversal 
indices behave trivially, so that the blocks we still have 
to diagonalise are quite small. The scalar block of P„ is 
given by 



P d 



-En, (79) 

0, (80) 

PnEn = a(r - dnn - dqq) + {1 - c)dqq, (81) 

P„r = {b-l)En, (82) 

the vector block by 



PndnA = —Ea, 
PndAn — 0) 
PuEa = —dnA + cdAn, 



(83) 
(84) 
(85) 



and the tensor block by PndAB = 0. The characteristic 
variables obey P„f7 = UA and comprise the four scalars 

dqq, (86) 
Uo = T+{b-l)dnn, (87) 
U± = a{T - dnn - dqq) + {1 - C)dqq ± VabEn, (SS) 



with speeds (0, 0, ±\/a6), six vector components 



(^An, 

U±A = dnA — cdAn T Ea, 



(89) 
(90) 



with speeds (0,±1), and three tensor components dAB 
with zero speed. The characteristic variables for the con- 
straints are 



Ca, C^ij , Cijk, 

C± = aCn + (1 - c)Cm^ ± V^Ce 



(91) 
(92) 



with speeds (0, ±va6). (We have not bothered to de- 
compose Cij and Cijk as all their components have zero 
speed.) Both the main and constraint system are strongly 
hyperbolic for ab > Q. Note that the coefficients of F in 
Uq etc. are themselves left (row) eigenvectors of P„.) The 
non-zero speed variables contain terms which are them- 
selves zero speed variables, (for example, U± contains 



dqq) 



In the base of characteristic variables the EM system 



is 



dtdqq = -d^{U+A-U-A), 

dtUo 



--{b-l)d^{U+A-U.A), 



±V^dnU± - ^^±^a^(L/. 



(93) 
(94) 

+A - U.a) 
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-Au+A 



+ U-a) + (c' - l)dAn 



(95) 



dtU±A = 



±dnU±A 



2Vab 
±d"{dBA -cdAs) 



dA{U+-U-) 



a6 + c — 1 



2ab 



(C/+ + C/_) + 



1 - c 



dtd-An 

dfdAB 



c-1 

' 2ab 
1 



2Vab 



(2c - 2 + 2a + ab)d, 
dA{U+ - C/_), 



99 



-79As5^((7+c-C/_c)- 



(96) 
(97) 

(98) 



Similarly, the constraint system in characteristic vari- 
ables is 

dtC± = ±Vd>dnC± 

±V^9^ [aCA + (1 - c)Ca/] , (99) 
Va6 , 



dtiaC'A) 
dtCij 
dtCijk 



-9^(C+-C_), 



0, 
0. 



(100) 

(101) 
(102) 



3. Symmetric hyperbolicity 

The first-order EM system admits the conserved en- 
ergy 

e = coeo+ciei, (103) 

where cq and ci are two arbitrary coefficients and 

eo = E'E, + d'^d,j ~cd'^dj,-2aTd,' 

+ {c-l-ab + 2a){d^)^, (104) 
£1 = [T+{b-l)dif. (105) 

The flux is 

F' ^ 2co [a(r - d/)£;* + (1 - c)d/£;* - {d'^ - cd^')Ej] 

(106) 

Note that ei has zero flux, e and i^" can be written in 
terms of characteristic variables as 

e = Co [dABd^"" - cdABd''^ + (1 - c2)d^„d-4„ 

+i(L/+^L/^ + L/_^c/:^) + ^(c/^ + ul) 

+ ^^(4a(igqC/o + (2 - 4a + a5 - 2c)d^,2) 



+ (ci-^co) [Uo + {b-l)dg 



(107) 



2^° 



(C/+AC/+ - U-aU^) 



1 



(108) 



Now we impose that e is positive definite. We can write 

^2 (109) 



dijd^-' — cdijd-''' ~ {cosipdij — sin(/3(iji) 



where sm2ip = c. This is positive definite for |c| < 1. 
Overall there are four positivity conditions: 



Co > 0, ci > 



3a' 



Sab + 2c-2 
2 

|c| < 1, ab > -(1 - c). 



Co, 



(110) 



The covariant constraint energy density and flux, after 
fixing an overall factor, are 

Cc - [aa + [1 - c)aj^]' + abCl 

+wi{ajf + W2{ajk)^ (111) 

= 2a6C£;[aC, + (l-c)C,/]. (112) 

It is unconditionally positive definite. In terms of char- 
acteristic variables, 



1 



pn 



-iCl + Cl) + [aCAHl-c)CA,^r 
\V^b{Cl-Cl). 



(113) 
(114) 



The terms multiplied by wi and W2 are strictly constant, 
and so wi > and W2 > are arbitrary. 

For the main system and the constraint system to be 
symmetric hyperbolic, the parameters of the system must 
obey the positivity conditions (|110|l . The strong hyper- 
bolicity condition ab > 0, which must also be obeyed, is 
implied by these. 



C. Pseudo-differential reduction to first order 

1. Field equations 

We now carry out the pseudo-differential reduction for 
comparison. Defining u as the Fourier transform of u in 
space, and substituting this into the field equations we 
obtain 

dtA, = -E,, (115) 
dtEi — ujjLO-' Ai + aiLOiV — (1 — a)ujiijj-' Aj^ (116) 
dtt - {b-l)iu,E\ (117) 

The fact that the system is second order in space is in- 
dicated by the fact that the highest power of \oj\ on the 
right-hand side is On the other hand, there is no 
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power of 1^1 on the right-hand side of dtAi. We define 
new variables by absorbing powers of into some of the 
Fourier transforms, namely 



(118) 



In these variables, the ODE system in the Fourier vari- 
ables becomes. 

dm = -t\uj\E„ (119) 

dfEi — i\Lu\ —di + aniT+(l — a)nin-'aj , (120) 

dtf = i\uj\{b - l)n,E\ (121) 

where again Ui = The right-hand sides are now 

the symbols of first-order pseudo-differential operators, 
but only the third is also the symbol of a differential 
operator. 

2. Strong hyperbolicity 

We can now find the characteristic variables of the 
pseudo-differential system using the same tensor decom- 
position as above. We write the pseudo-differential sys- 
tem as 



dtii = iliulPnU. 



(122) 



P„ in (|122|l is a submatrix of P„ in lf75|l . with Ei ^ Ei, 
r r, dnAi — > = di, but dA^i because there 

are no transversal derivatives in the pseudo-differential 
framework. Therefore the rows and columns correspond- 
ing to dAi and dqq disappear. Note that with these dele- 
tions, c disappears from P„: this must happen as c pa- 
rameterises an ambiguity of the differential reduction. 
The characteristic variables are 



Uo = r+(6-l)a„, 
tj± = a(f — d) ± VabEn, 
U±B = cibTEb, 



(123) 
(124) 
(125) 



with speeds (0,±vafe, ±1). 

The Fourier transforms of the constraints are 



Ce - Cr=t- a„. (126) 

With 

cr = i\uj\Cr (127) 

we have 

dtCs = i\uj\acr, (128) 

dtcr = i\uj\bCE. (129) 

The characteristic variables are 

C± = acr±V^CE^i\uj\U± (130) 



with speeds iyab. The pseudo-differential system is 
strongly hyperbolic for ah > 0. 

Because of the relation of the two matrices P„ dis- 
cussed above, the characteristic variables for the pseudo- 
differential reduction could be obtained from those for 
the differential reduction by the replacements dAi 0, 
dni di, Ei Ei, r ^ F, and Ce ^ Ce, Ca 0, 



Cn > Cr, Ci 



0, C, 



ijk 



0. 



D. Second-order in space, first order in time 

1. Strong hyperbolicity 

To construct the second-order characteristic variables, 
we consider the principal part of the evolution of u = 
{Ei ,diAj,r), which obey 



dtU = PndnU + P OaU. 



(131) 



One might think that P„ is the same matrix as in the dif- 
ferential first order reduction. But consider, for example, 
Eq. H85|l . which is shorthand for 

dtEA — —dndnA + cdndAn + transversal derivatives 

(132) 

In the second-order system this becomes 

dtEA — —O^Aa + cdndAAn + transversal derivatives, 

(133) 

but we can commute derivatives and so lump dndAAn = 
dAdnAn with the transversal derivatives. As dtdAAn = 
dAdtAn can also be considered purely transversal, we can 
set both the row and column corresponding to OaA^ in 
P„ equal to zero. This is true for all OaAi, and dACr in 
the constraint system. Therefore we are effectively diag- 
onalising smaller matrices P„ for the main and constraint 
systems. These are identical (up to variable names) to 
those in the pseudo-differential reduction. 

The characteristic variables for the second-order EM 
system obtained by diagonalising this reduced P„ are 



Uo = T + (b - l)dr,An, 
U'± = a{T - dnAn) ± Vd>Er,, 



(134) 
(135) 
(136) 

with speeds (0, ±\/a6, ±1) for the main system, and 

C± = ad„Cr ± VabCE (137) 



with speeds iyab for the constraint system. To complete 
the set of characteristic variables, we add the transversal 
derivatives 

OaA, dACr, (138) 
which are by definition zero-speed variables. 
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2. Symmetric hyperbolicity 

We obtain a conserved energy for the second-order sys- 
tem by writing down the most general scalar e that is 
quadratic in F, Ei and diAj and adjusting the free coef- 
ficients to obtain conservation. The result is 



ciei + €262, 

2 - 2aTd^A, 



e = coeo 

+ {2a-l~ab){d'A,f 
61 = [r+(6-l)aM,]2, 
£2 - id'A,f-d,Ajd^A\ 



(139) 

(140) 
(141) 
(142) 



The flux is given by 
= 2{[(a7 



- ^xiy^,^il~a)d^A,]E^ 
= 2{d'A'Ej-d'AjE'), 



d'A^Ej} , (143) 
(144) 



while ei has no flux. The constraint energy is unique up 
to an overall factor, and is given by 

= a'^d^Crd'Cr + abCl, (145) 
F^ = 2a%Crd'CE. (146) 

If we add zero-speed terms OaAi to Uj. and U'j_j^ as 
follows, 



U± 

U±B 



dnA„ - d^As) ± VabE, 



a(r^ 

-1-(1-c)9^Ab, 
dnAs T Eb ~ cdsAr, 



(147) 
(148) 



where c = C2/C0, we can write the main flux in the form 
(|l()8|l . and with (|137|l we can write the constraint flux 
in the form (|114|l . Furthermore, the evolution equations 
expressed in the second-order characteristic variables are 
given by H93I98|I . with the obvious replacements. This 
works because the same constant c arises twice in differ- 
ent ways: in the flrst-order system it parameterises an 
arbitrariness of writing the second-order field equations 
in terms of first-order variables, and this is passed on to 
the energy. In the second-order system there is no such 
ambiguity in the field equations, but c arises as the free 
coeflScient C2/C0 in the energy. 



E. Constraint-preserving boundary conditions 

The following discussion applies equally to the flrst- 
order and second-order symmetric hyperbolic systems. 
The constraint-preserving boundary conditions we are 
about to derive contain only time derivatives when ap- 
plied to the first-order system, but both space and time 
derivatives when applied to the second-order system. We 
shall explicitly write down expressions for the first-order 
system. The corresponding expressions for the second or- 



der system are obtained by the replacement d^j — > diA 
and interpreting the characteristic variables as those of 
the second-order system. 



The maximally dissipative boundary conditions on the 
main variables are 



U+A - K2U-A = Ia, 



(149) 
(150) 



where |ki| < 1 and |k2| < 1, and / and Ja are given func- 
tions. We are only interested in solutions with vanishing 
constraints. Therefore we formally impose homogeneous 
maximally dissipative boundary conditions on the con- 
straint system, that is 



K3C- = 



(151) 



with IK3I < 1. We now show how these boundary condi- 
tions can be made consistent, following the general pro- 
cedure set out in Sec. [V| 

C± can be expressed in terms of the main characteristic 
variables as 

1 



tV^{U+a - U-a) 

+ {c-l){U+A + U-A + 2cdAn)]. (152) 



We use the evolution equation for U± to replace dnU± 
by dtU±. Eq. (|151|l then becomes 



dt{U++K3U-) 

(1 - K3)\fab 



(153) 



d''[U+A + U-A + 2{c-l)dAn]- 

We impose this as an evolution equation for an auxiliary 
field X = U+ + on the boundary. We then im- 

pose the boundary condition p49|l on C/-|_ with ki = — K3 
and f = X. We also impose the maximally dissipative 
boundary conditions 1)1501) on U+a, which represent the 
two polarisations of ingoing radiation. K2 and /a must 
therefore be chosen on physical grounds: stability only 
dictates |k2| < 1. 



1. Dirichlet boundary conditions 

As discussed in general in Sec.|3 one can demonstrate 
well-posedness with these constraint-preserving bound- 
ary conditions if the bulk system does not couple back 
to the main system. There are two cases in which this 
can be achieved: For K3 = 1, dAn is eliminated from the 
boundary system, while for K3 = —1, its evolution equa- 
tion H97(l becomes part of an expanded closed boundary 
system. For K3 = K2 = —ki = 1, the boundary system 
and boundary conditions are 

dtX = (1-c)9^/a, (154) 
U+ + U+ = X, (155) 
U+A-U-A = Ja, (156) 

where /a is free data. The boundary condition on U± 
corresponds to Dirichlet boundary conditions on C^;. 
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2. Neumann boundary conditions 

For K3 = K2 = — = —1, the boundary system and 
the boundary conditions for the main system are 

dtX = ^d^[fA + 2{c-l)dAn], (157) 



dtdAn = r^dAX, (158) 

2vao 

U+-U- = X, (159) 

U+A + U-A - fA. (160) 



The boundary condition on U± corresponds to Neumann 
boundary conditions on Cr- 

The boundary system is symmetric hyperbohc with the 
energy density, flux, and source term 

eb = X^ +iab{l-c){dAnf, (161) 

dtEb = dAF^ + Sb, (162) 

= 4\/^(c- l)XdA„, (163) 

Sb = 2VabXd'^fA. (164) 

The characteristic variables of the boundary system in 
the direction are 

X± =XT2Va^(l-c)d™„, (165) 

with speeds ±\/l — c. The boundary energy is controUed 

by 

^ = / sdS+ f F^ds, (166) 

Fr = IVT^^Xl-Xl), (167) 

where myi is the normal to ddfl in dfl. 

For a smooth boundary {ddfl — 0), the boundary en- 
ergy can therefore be estimated by 

^ < 2V^b\\X\\ \\d^fA\\ < 2V^by/E'b\\d^fA\\ (168) 
dt 

For a cubic boundary, one also has to consider the ddfl 
terms at the edges of each face. From (|160|) and (|90|) we 
have 

dmn = ^ 5 (169 

1 — 

where /j^"'' are the free data on the boundary with normal 
Ui. They are already given. This leaves us no choice but 
to impose Dirichlet boundary conditions given by 1)169(1 
on d„in at the edge with normal niA of the face with nor- 
mal rii. From (|165|l and 1)1671) we see that these are max- 
imally dissipative boundary conditions on the system on 
each face, so that the estimate H168|) still holds. There- 
fore the coupled bulk/boundary system is well-posed in 
the Neumann case. 



3. Laplace- Fourier analysis 

With our current methods we can decouple the bound- 
ary system from the bulk system, and hence prove well- 
posedness of the initial-boundary value problem, only for 
the two extreme values K3 = ±1. The same was found 
previously in imposing constraint-preserving boundary 
conditions on the linearised Einstein-Christoffel system 
I3 and on the full Einstein equations in harmonic gauge 
^ . But the lack of a proof does not rule out that for ex- 
ample the Sommerfeld boundary conditions with K3 = 
could also be well-posed. 

To investigate this possibility further, we now use the 
Laplace-Fourier method to explicitly check for modes 
that grow arbitrarily rapidly. The EM equations are 
untypical in that they are linear with constant coeffi- 
cients. For a quasilinear system such as the Einstein 
equations, we use the high frequency, or frozen coeffi- 
cients, approximation, and consider only the principal 
part of the evolution equations and constraints, thus re- 
covering a linear system with constant coefficients. In 
the frozen coefficients approximation we assume that the 
linearised perturbation varies over space and time scales 
much smaller than those given by the background solu- 
tion and the numerical domain. For consistency we must 
therefore assume that the domain is a half space and that 
the boundary is a plane. Without loss of generality, we 
assume that the evolution domain is — 00 < < and 
— 00 < = (a;^,x^) < 00. The general solution with 
homogeneous boundary conditions can then be written 
as a sum of modes of the form 

u{x\ t) = e^t+^^^^-'+^^'i; (170) 

where w is a constant eigenvector to be determined. (If 
the corresponding eigenvalue A was degenerate, v would 
be a polynomial in rather than constant, but this does 
not happen in for the EM system.) If the boundary con- 
ditions are inhomogeneous, we must add to this ansatz a 
particular integral that does not concern us here because 
its growth is controlled by the boundary data. We are 
interested only in modes that are square-integrable over 
space at any moment in time. Therefore we assume that 
LOA is real, and that ReA > 0. s and A will in general be 
complex. 

If a mode of this form exists for some (s, uja^ A, v), then 
one exists also for (fcs, kojAi kX, v) for any fc > 0. There- 
fore, if any growing mode exists, there are growing modes 
with arbitrarily large growth rates ks, and the problem 
is ill-posed. A necessary condition for well-posedness of 
the initial-boundary value problem is therefore that the 
homogeneous boundary conditions rule out the existence 
of any mode with Res > for all real uja and all A with 
ReA > 0. 

We now examine our proposed constraint-preserving 
boundary conditions for the EM system. To simplify the 
algebra, we restrict ourselves to the case a = b = 1, 
c = (which is symmetric hyperbolic). The evolution 
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equations for the characteristic variables simplify to 



dfdqq 

dtUo 

dtdAn 

dfdAB 
dtU±A 



d^{U+A - U^a), 



--dA{U+ -[/_), 



(171) 
(172) 
(173) 



±dnU±Td'^dAn, 



(175) 



±dnU±A ± d^'dBA T OaUo ± T;d^dqq. (176) 



(Once again, we could trivially rewrite this for the 
second-order system) . The constraint-preserving bound- 
ary conditions are 



+A 



Kali- A 



~il-K3)dAn\, (177) 

U+A-H2U-A = 0. (178) 

Using the Laplace-Fourier ansatz in the evolution equa- 
tions, we obtain a linear eigenvalue problem for v with 
eigenvalue A(s, \u!\). Its solution is 



[/_ 

Ua+ 
Ua- 



7ie' 



(179) 



where A — -^/s^ + Iwp and if) — yzM-T — z with z = 
s/\ijo\. 7i and 72A are free coefhcients. We ensure that 
ReA > for Res > by putting the branch cut for the 
square root on the negative real axis. 

Applying the boundary condition (|178f) to 118UI) gives 
the algebraic condition 72^1 ('«2'0^ — 1) = 0. This means 
that 72A = unless K2i'^{z) — 1. Now i^^lz) maps the 
right half-plane into the unit disk with the negative real 
axis removed. Therefore, k,2'4'^{z) = 1 for Rez > is 
ruled out if |k2| < 1. The boundary condition 1)1 77(1 
applied to H18U|) gives a more complicated algebraic con- 
dition linking 71 and u)^"f2A- However, with 72A = this 
simplifies, after some algebra, to 71(^3?/'^ ^ 1) = 0, and 
so we must also impose l^al < 1. 

Our calculation has shown that the initial-boundary 
value problem has no growing square-integrable modes 
of the form ^17i)ll80|l for |k2| < 1 and l^tgl < 1. This 
is consistent with our energy method proof that the EM 
initial-boundary value problem is well-posed for K3 = ± 1 
and |k2| < 1, but it also suggests that the problem is 
well-posed in the more general case, and we expect that 
a proof of this can be found using a modified energy. 



VII. CONCLUSIONS 

We have proposed a new method for setting up a well- 
posed initial-boundary value problem for evolution sys- 



tems, such as variants of the ADM equations, that are 
first order in time but second order in space. One might 
call it the second-order energy method. It involves two 
steps: 

1. Find a complete set of characteristic variables for 
any given direction that obey Eq. H34|l . They are not yet 
unique. This is our definition of strong hyperbolicity. 

2. Find a conserved covariant energy, and express it 
and the corresponding flux in terms of the characteristic 
variables. This can be done for a unique set of char- 
acteristic variables. This is our definition of symmetric 
hyperbolicity. 

Our definition of strong hyperbolicity for a second- 
order system is equivalent to that using a pseudo- 
differential reduction to first order P, Q, as both def- 
initions require the same matrix to be diagonalisable. 
But the pseudo-differential method cannot be extended 
to symmetric hyperbolicity, roughly speaking because 
Fourier transforms do not allow for boundaries in space. 
Our definition of symmetric hyperbolicity for a second- 
order system is similar to that for the differential reduc- 
tion to first order, but it does not introduce auxiliary 
variables and constraints, and so does not enlarge the 
solution space. 

We have worked through both types of first-order re- 
duction, as well as the second-order energy method, in 
the example of electromagnetism, in order to illustrate 
the rather subtle similarities and differences. In a com- 
panion paper [l^ we shall use the energy method to 
prove symmetric hyperbolicity and suggest constraint- 
preserving boundary conditions for variants of the ADM 
equations. 
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APPENDIX A: AN ALTERNATIVE 
FIRST-ORDER REDUCTION 

Lindblom et al |ll| have imposed constraint-preserving 
boundary conditions of the Sommerfeld type on a sym- 
metric hyperbolic first-order reduction of the EM system 
(without proving well-posedness of the resulting initial- 
boundary value problem). They consider the following 
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first-order reduction of the EM system: 

dtE, = - rf,,) + ^{dd/ - a^dy), (Ai) 

dtd,, = -^,E,+J2S^J^''Ek. (A2) 

This corresponds to our first-order reduction with a — 

(and r therefore not part of the system) and c — 

1 — (71/2), plus the additional term parameterised by 
72. The latter plays a similar role to our parameter 
6, in adding the constraint Cb to the evolution equa- 
tion for di'i— r). They find that the system is strongly 



hyperbolic for 7172 > 0, and symmetric hyperbolic for 
< 7i < 4 and 72 > 1/3. However, this first-order reduc- 
tion has as its second-order counterpart only the original 
EM system (|6UI61ll . which is only weakly hyperbolic. The 
term multiplied by 71 vanishes identically because par- 
tial derivatives can be commuted, and the term 72 cannot 
arise if we evolve Ai rather than dij. In hindsight, one 
can see that both terms have second-order counterparts 
as long as we have the trace of dij as an auxiliary variable: 
the rest of is not required. 
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